/*
Copyright 2013 Lin, Chung-Ming (aka Alva Lin. E-mail: alva0930@gmail.com)

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/
import java.io.*;
import java.util.*;
import Jama.*;	// The 'Jama' package could be downloaded from: http://math.nist.gov/javanumerics/jama/
import org.apache.hadoop.conf.Configuration;
import org.apache.hadoop.fs.*;
import org.apache.hadoop.io.*;
import org.apache.hadoop.io.compress.*;

/**
 * A collection of utility functions for {@link Matrix} class.
 *
 * @see	Matrix
 * @author	Lin, Chung-Ming (aka Alva Lin. E-mail: alva0930@gmail.com)
 */
public class MatrixFunction
{
	private static final double EPS = 1e-12;

/* ------------------------
   Public Methods
 * ------------------------ */

	/**
	 * Horizontally append {@link Matrix} B to {@link Matrix} A.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}} and {@link Matrix} B represents {{5.0, 6.0}, {7.0, 8.0}},
	 * then {@link Matrix} C = appendH(A, B) represents {{1.0, 2.0, 5.0, 6.0}, {3.0, 4.0, 7.0, 8.0}}.
	 * The number of rows of {@link Matrix} A must be the same as the number of rows of {@link Matrix} B,
	 * otherwise it will cause an {@link AssertionError} when assertions are enabled,
	 * or cause unexpected result when assertions are disabled.
	 *
	 * @param	A	the {@link Matrix} to be put on the left
	 * @param	B	the {@link Matrix} to be put on the right
	 * @return	the {@link Matrix} formed by horizontally append {@link Matrix} B to {@link Matrix} A
	 */
	public static Matrix appendH(Matrix A, Matrix B)
	{
		int nRow1 = A.getRowDimension();
		int nCol1 = A.getColumnDimension();
		int nRow2 = B.getRowDimension();
		int nCol2 = B.getColumnDimension();

		assert(nRow1 == nRow2);

		Matrix output = new Matrix(nRow1, nCol1 + nCol2);
		output.setMatrix(0, nRow1 - 1, 0, nCol1 - 1, A);
		output.setMatrix(0, nRow1 - 1, nCol1, nCol1 + nCol2 - 1, B);

		return output;
	}

	/**
	 * Read a {@link Matrix} from a csv(comma-separated values) file.
	 *
	 * @param	path	the path of the csv file
	 * @return	the {@link Matrix} read from the csv file
	 */
	public static Matrix csvread(String path) throws FileNotFoundException, IOException
	{
		return csvread(path, false);
	}

	/**
	 * Read a {@link Matrix} from a csv(comma-separated values) file.
	 *
	 * @param	path	the path of the csv file
	 * @param	skipFirstLine	whether or not to skip the first line of the csv file
	 * @return	the {@link Matrix} read from the csv file
	 */
	public static Matrix csvread(String path, boolean skipFirstLine) throws FileNotFoundException, IOException
	{
		BufferedReader input = new BufferedReader(new FileReader(path));

		String record = input.readLine();
		if (skipFirstLine) {record = input.readLine();}
		double[] numbers = getDoubleArray(record);
		ArrayList<Double>[] lists = new ArrayList[numbers.length];
		for (int i = 0; i < lists.length; i++)
		{
			lists[i] = new ArrayList<Double>();
			lists[i].add(numbers[i]);
		}

		while ((record = input.readLine()) != null)
		{
			numbers = getDoubleArray(record);
			for (int i = 0; i < lists.length; i++)
			{
				lists[i].add(numbers[i]);
			}
		}

		input.close();

		return arrayListsToMatrixTranspose(lists);
	}

	/**
	 * Write a {@link Matrix} to a csv(comma-separated values) file.
	 *
	 * @param	path	the path of the csv file
	 * @param	A	the {@link Matrix} to be written
	 */
	public static void csvwrite(String path, Matrix A) throws IOException
	{
		csvwrite(path, A, false);
	}

	/**
	 * Write a {@link Matrix} to a csv(comma-separated values) file.
	 *
	 * @param	path	the path of the csv file
	 * @param	A	the {@link Matrix} to be written
	 * @param	convertToInt	whether or not to convert the contents of {@link Matrix} A to integer
	 */
	public static void csvwrite(String path, Matrix A, boolean convertToInt) throws IOException
	{
		BufferedWriter output = new BufferedWriter(new FileWriter(path));

		for (double[] da : A.getArray())
		{
			StringBuffer sb = null;
			for (double d : da)
			{
				if (null == sb)
				{
					sb = new StringBuffer();
				}
				else
				{
					sb.append(",");
				}

				if (convertToInt)
				{
					sb.append("" + ((int)d));
				}
				else
				{
					sb.append("" + d);
				}
			}
			output.write(sb.toString());
			output.newLine();
		}

		output.close();
	}

	/**
	 * Compare each element of a {@link Matrix} to a scalar value.
	 * 1 indicates the element in the corresponding position is the same as the scalar value, 0 indicates they are different.
	 * For example, assume that {@link Matrix} A represents {{2.0, 3.0}, {3.0, 4.0}} and the value is 3.0,
	 * then {@link Matrix} B = equal(A, value) represents {{0.0, 1.0}, {1.0, 0.0}}.
	 *
	 * @param	A	a {@link Matrix} to be compared with the scalar value
	 * @param	value	a scalar value to be compared with the {@link Matrix}
	 * @return	the result {@link Matrix} of comparison between A and value
	 */
	public static Matrix equal(Matrix A, double value)
	{
		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();

		Matrix output = new Matrix(nRow, nCol);
		for (int i = 0; i < nRow; i++)
		{
			for (int j = 0; j < nCol; j++)
			{
				if (Math.abs(A.get(i, j) - value) < EPS)
				{
					output.set(i, j, 1);
				}
			}
		}

		return output;
	}

	/**
	 * Element-by-element comparison between two {@link Matrix} objects.
	 * 1 indicates the two elements in the same position are the same, 0 indicates they are different.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}} and {@link Matrix} B represents {{1.0, 6.0}, {7.0, 4.0}},
	 * then {@link Matrix} C = equal(A, B) represents {{1.0, 0.0}, {0.0, 1.0}}.
	 * The number of rows of {@link Matrix} A must be the same as the number of rows of {@link Matrix} B,
	 * and the number of columns of {@link Matrix} A must be the same as the number of columns of {@link Matrix} B,
	 * otherwise it will cause an {@link AssertionError} when assertions are enabled,
	 * or cause unexpected result when assertions are disabled.
	 *
	 * @param	A	a {@link Matrix} to be compared with the other {@link Matrix}
	 * @param	B	a {@link Matrix} to be compared with the other {@link Matrix}
	 * @return	the result {@link Matrix} of the element-by-element comparison between A and B
	 */
	public static Matrix equal(Matrix A, Matrix B)
	{
		int nRow1 = A.getRowDimension();
		int nCol1 = A.getColumnDimension();
		int nRow2 = B.getRowDimension();
		int nCol2 = B.getColumnDimension();

		assert((nRow1 == nRow2) && (nCol1 == nCol2));

		Matrix output = new Matrix(nRow1, nCol1);
		for (int i = 0; i < nRow1; i++)
		{
			for (int j = 0; j < nCol1; j++)
			{
				if (Math.abs(A.get(i, j) - B.get(i, j)) < EPS)
				{
					output.set(i, j, 1);
				}
			}
		}

		return output;
	}

	/**
	 * Read a {@link Matrix} from a <a href="http://wiki.apache.org/hadoop/SequenceFile">sequence file</a> written by hdfswrite(...) on HDFS(Hadoop Distributed File System).
	 *
	 * @param	path	the path of the file on HDFS
	 * @return	the {@link Matrix} read from the file
	 */
	public static Matrix hdfsread(Path path) throws IOException
	{
		Configuration conf = new Configuration();
		FileSystem fs = FileSystem.get(conf);
		IntWritable key = new IntWritable();
		Text value = new Text();
		SequenceFile.Reader input = new SequenceFile.Reader(fs, path, conf);

		input.next(key, value);
		double[] numbers = getDoubleArray(value.toString());
		ArrayList<Double>[] lists = new ArrayList[numbers.length];
		for (int i = 0; i < lists.length; i++)
		{
			lists[i] = new ArrayList<Double>();
			lists[i].add(numbers[i]);
		}

		while (input.next(key, value))
		{
			numbers = getDoubleArray(value.toString());
			for (int i = 0; i < lists.length; i++)
			{
				lists[i].add(numbers[i]);
			}
		}

		input.close();

		return arrayListsToMatrixTranspose(lists);
	}

	/**
	 * Write a {@link Matrix} to a <a href="http://wiki.apache.org/hadoop/SequenceFile">sequence file</a> on HDFS(Hadoop Distributed File System).
	 * Each row is written as a record.
	 * The key({@link IntWritable}) of the record is the index (start from 0) of the row.
	 * The value({@link Text}) of the record is the text of the comma-separated values of the row.
	 *
	 * @param	path	the path of the file on HDFS
	 * @param	A	the {@link Matrix} to be written
	 */
	public static void hdfswrite(Path path, Matrix A) throws IOException
	{
		hdfswrite(path, A, SequenceFile.CompressionType.NONE, null);
	}

	/**
	 * Write a {@link Matrix} to a <a href="http://wiki.apache.org/hadoop/SequenceFile">sequence file</a> on HDFS(Hadoop Distributed File System).
	 * Each row is written as a record.
	 * The key({@link IntWritable}) of the record is the index (start from 0) of the row.
	 * The value({@link Text}) of the record is the text of the comma-separated values of the row.
	 *
	 * @param	path	the path of the file on HDFS
	 * @param	A	the {@link Matrix} to be written
	 * @param	type	the compression type
	 * @param	codec	the compression codec
	 */
	public static void hdfswrite(Path path, Matrix A, SequenceFile.CompressionType type, CompressionCodec codec) throws IOException
	{
		Configuration conf = new Configuration();
		FileSystem fs = FileSystem.get(conf);
		IntWritable key = new IntWritable();
		Text value = new Text();
		SequenceFile.Writer output = SequenceFile.createWriter(fs, conf, path, key.getClass(), value.getClass(), type, codec);

		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();
		double[][] matrix = A.getArray();

		for (int i = 0; i < nRow; i++)
		{
			StringBuffer sb = new StringBuffer("" + matrix[i][0]);
			for (int j = 1; j < nCol; j++)
			{
				sb.append("," + matrix[i][j]);
			}
			key.set(i);
			value.set(sb.toString());
			output.append(key, value);
		}

		output.close();
	}

	/**
	 * Get the log value of each element of a {@link Matrix}.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}},
	 * then {@link Matrix} B = log(A) represents {{{@link Math}.log(1.0), {@link Math}.log(2.0)}, {{@link Math}.log(3.0), {@link Math}.log(4.0)}}.
	 *
	 * @param	A	a {@link Matrix} to apply the log function to
	 * @return	the result {@link Matrix} of applying the log function to each element of A
	 */
	public static Matrix log(Matrix A)
	{
		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();

		Matrix output = new Matrix(nRow, nCol);
		for (int i = 0; i < nRow; i++)
		{
			for (int j = 0; j < nCol; j++)
			{
				output.set(i, j, Math.log(A.get(i, j)));
			}
		}

		return output;
	}

	/**
	 * Get the index of the (first) maximum element of each row.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0, 3.0}, {3.0, 4.0, 4.0}, {5.0, 4.0, 3.0}},
	 * then {@link Matrix} B = maxIndexRow(A) represents {{2.0}, {1.0}, {0.0}}.
	 *
	 * @param	A	a {@link Matrix} to be looked for the index of the (first) maximum element of each row
	 * @return	a single column {@link Matrix} composed of the index of the (first) maximum element of each row in {@link Matrix} A
	 */
	public static Matrix maxIndexRow(Matrix A)
	{
		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();

		Matrix output = new Matrix(nRow, 1);
		for (int i = 0; i < nRow; i++)
		{
			int maxIndex = 0;
			for (int j = 1; j < nCol; j++)
			{
				if (A.get(i, j) > (A.get(i, maxIndex) + EPS))
				{
					maxIndex = j;
				}
			}
			output.set(i, 0, maxIndex);
		}

		return output;
	}

	/**
	 * Get the mean of elements of a {@link Matrix}.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}},
	 * then double m = mean(A) represents 2.5.
	 *
	 * @param	A	a {@link Matrix} to be calculated for the mean
	 * @return	the mean of elements of {@link Matrix} A
	 */
	public static double mean(Matrix A)
	{
		return (sum(A) / (A.getRowDimension() * A.getColumnDimension()));
	}

	/**
	 * Subtract each element of a {@link Matrix} by a scalar value.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}},
	 * then {@link Matrix} B = minus(A, 2.0) represents {{-1.0, 0.0}, {1.0, 2.0}}.
	 *
	 * @param	A	a {@link Matrix} to be subtracted by the scalar value
	 * @param	value	a scalar value used as subtrahend
	 * @return	the result {@link Matrix} of each element of A subtracted by value
	 */
	public static Matrix minus(Matrix A, double value)
	{
		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();

		Matrix B = new Matrix(nRow, nCol, value);
		
		return A.minus(B);
	}

	/**
	 * Subtract a scalar value by each element of a {@link Matrix}.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}},
	 * then {@link Matrix} B = minus(2.0, A) represents {{1.0, 0.0}, {-1.0, -2.0}}.
	 *
	 * @param	A	a {@link Matrix} used as subtrahend
	 * @param	value	a scalar value to be subtracted by each element of the {@link Matrix}
	 * @return	the result {@link Matrix} of value subtracted by each element of A
	 */
	public static Matrix minus(double value, Matrix A)
	{
		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();

		Matrix B = new Matrix(nRow, nCol, value);
		
		return B.minus(A);
	}

	/**
	 * Predict the classes (ranged from 0 to &ltnumber of labels - 1&gt) of input features by a single-hidden-layer logistic-units neural network.
	 * The number of columns of {@link Matrix} theta1 must be equal to the number of columns of {@link Matrix} X plus 1,
	 * and the number of columns of {@link Matrix} theta2 must be equal to the number of rows of {@link Matrix} theta1 plus 1,
	 * otherwise it will cause an {@link AssertionError} when assertions are enabled,
	 * or cause unexpected result when assertions are disabled.
	 *
	 * @param	X	the features of the test cases, of size &ltnumber of test cases&gt X &ltnumber of features (input layer size)&gt
	 * @param	theta1	weights for the input layer to the hidden layer, of size &lthidden layer size&gt X &ltinput layer size + 1&gt
	 * @param	theta2	weights for the hidden layer to the output layer, of size &ltnumber of labels (output layer size)&gt X &lthidden layer size + 1&gt
	 * @return	the {@link Matrix} of the classification prediction, of size &ltnumber of test cases&gt X 1
	 */
	public static Matrix predictSingleHiddenLayerNN(Matrix X, Matrix theta1, Matrix theta2)
	{
		assert((X.getColumnDimension() + 1) == theta1.getColumnDimension());
		assert((theta1.getRowDimension() + 1) == theta2.getColumnDimension());
		int nRow = X.getRowDimension();

		Matrix h1 = sigmoid(appendH(new Matrix(nRow, 1, 1), X).times(theta1.transpose()));
		Matrix h2 = sigmoid(appendH(new Matrix(nRow, 1, 1), h1).times(theta2.transpose()));

		return maxIndexRow(h2);
	}

	/**
	 * Get a {@link Matrix} with uniformly distributed random elements.
	 *
	 * @param	nRow	number of rows of the Matrix, must be greater than 0
	 * @param	nCol	number of columns of the Matrix, must be greater than 0
	 * @param	range	the generated random elements will be between range and (-range)
	 * @return	a {@link Matrix} with uniformly distributed random elements between range and (-range)
	 */
	public static Matrix random(int nRow, int nCol, double range)
	{
		assert(nRow > 0);
		assert(nCol > 0);

		Random r = new Random();
		Matrix output = new Matrix(nRow, nCol);
		for (int i = 0; i < nRow; i++)
		{
			for (int j = 0; j < nCol; j++)
			{
				output.set(i, j, r.nextDouble() * 2 * range - range);
			}
		}

		return output;
	}

	/**
	 * Reshape the elements of a range of an array to a {@link Matrix}.
	 * For example, assume that an array values represents {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0},
	 * then {@link Matrix} A = reshape(values, 0, 3, 2, 2) represents {{1.0, 3.0}, {2.0, 4.0}},
	 * and {@link Matrix} B = reshape(values, 4, 12, 3, 3) represents {{5.0, 8.0, 11.0}, {6.0, 9.0, 12.0}, {7.0, 10.0, 13.0}}.
	 * The length of the range (end - begin + 1) must be the same as the size of the {@link Matrix} (nRow * nCol),
	 * otherwise it will cause an {@link AssertionError} when assertions are enabled,
	 * or cause unexpected result when assertions are disabled.
	 *
	 * @param	values	an array to be reshaped to a {@link Matrix}
	 * @param	begin	the begin index of the array to be reshaped
	 * @param	end	the end index of the array to be reshaped
	 * @param	nRow	number of rows of the output {@link Matrix}, must be greater than 0
	 * @param	nCol	number of columns of the output {@link Matrix}, must be greater than 0
	 * @return	the {@link Matrix} reshaped from a range of an array
	 */
	public static Matrix reshape(double[] values, int begin, int end, int nRow, int nCol)
	{
		assert(nRow > 0);
		assert(nCol > 0);
		assert((end - begin + 1) == (nRow * nCol));

		Matrix output = new Matrix(nRow, nCol);
		for (int i = 0; i < nCol; i++)
		{
			for (int j = 0; j < nRow; j++)
			{
				output.set(j, i, values[begin + i * nRow + j]);
			}
		}

		return output;
	}

	/**
	 * Get the sigmoid value of each element of a {@link Matrix}.
	 * For example, assume that S(t) = (1 / (1 + {@link Math}.exp(-t))) and {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}},
	 * then {@link Matrix} B = sigmoid(A) represents {{S(1.0), S(2.0)}, {S(3.0), S(4.0)}}.
	 *
	 * @param	A	a {@link Matrix} to apply the sigmoid function to
	 * @return	the result {@link Matrix} of applying the sigmoid function to each element of A
	 */
	public static Matrix sigmoid(Matrix A)
	{
		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();

		Matrix output = new Matrix(nRow, nCol);
		for (int i = 0; i < nRow; i++)
		{
			for (int j = 0; j < nCol; j++)
			{
				output.set(i, j, sigmoid(A.get(i, j)));
			}
		}

		return output;
	}

	/**
	 * Get the sigmoid gradient value of each element of a {@link Matrix}.
	 * For example, assume that S(t) = (1 / (1 + {@link Math}.exp(-t))), SG(t) = (S(t) * (1 - S(t)))
	 * and {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}},
	 * then {@link Matrix} B = sigmoidGradient(A) represents {{SG(1.0), SG(2.0)}, {SG(3.0), SG(4.0)}}.
	 *
	 * @param	A	a {@link Matrix} to apply the sigmoid gradient function to
	 * @return	the result {@link Matrix} of applying the sigmoid gradient function to each element of A
	 */
	public static Matrix sigmoidGradient(Matrix A)
	{
		Matrix sA = sigmoid(A);
		return sA.arrayTimes(minus(1, sA));
	}

	/**
	 * Get the sum of elements of a {@link Matrix}.
	 * For example, assume that {@link Matrix} A represents {{1.0, 2.0}, {3.0, 4.0}},
	 * then double s = sum(A) represents 10.0.
	 *
	 * @param	A	a {@link Matrix} to be calculated for the sum
	 * @return	the sum of elements of {@link Matrix} A
	 */
	public static double sum(Matrix A)
	{
		int nRow = A.getRowDimension();
		int nCol = A.getColumnDimension();

		double output = 0;
		for (int i = 0; i < nRow; i++)
		{
			for (int j = 0; j < nCol; j++)
			{
				output += A.get(i, j);
			}
		}

		return output;
	}

	/**
	 * Unroll the elements of two {@link Matrix} objects to a array.
	 * For example, assume that {@link Matrix} A represents {{1.0, 3.0}, {2.0, 4.0}},
	 * and {@link Matrix} B represents {{5.0, 8.0, 11.0}, {6.0, 9.0, 12.0}, {7.0, 10.0, 13.0}},
	 * then double[] da = unroll(A, B) represents {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0}.
	 *
	 * @param	A	a {@link Matrix} to be unrolled in the left part of output array
	 * @param	B	a {@link Matrix} to be unrolled in the right part of output array
	 * @return	the array of elements of {@link Matrix} A unrolled in the left and elements of {@link Matrix} B unrolled in the right
	 */
	public static double[] unroll(Matrix A, Matrix B)
	{
		int nRow1 = A.getRowDimension();
		int nCol1 = A.getColumnDimension();
		int size1 = nRow1 * nCol1;
		int nRow2 = B.getRowDimension();
		int nCol2 = B.getColumnDimension();
		int size2 = nRow2 * nCol2;

		double[] output = new double[size1 + size2];
		for (int i = 0; i < nCol1; i++)
		{
			for (int j = 0; j < nRow1; j++)
			{
				output[i * nRow1 + j] = A.get(j, i);
			}
		}
		for (int i = 0; i < nCol2; i++)
		{
			for (int j = 0; j < nRow2; j++)
			{
				output[size1 + i * nRow2 + j] = B.get(j, i);
			}
		}

		return output;
	}

	private static Matrix arrayListsToMatrixTranspose(ArrayList<Double>[] lists)
	{
		int R = lists[0].size();
		int C = lists.length;
		double[][] matrix = new double[R][C];
		for (int i = 0; i < R; i++)
		{
			for (int j = 0; j < C; j++)
			{
				matrix[i][j] = lists[j].get(i);
			}
		}
		return new Matrix(matrix);
	}

	private static double[] getDoubleArray(String line)
	{
		String[] strings = line.split(",", 0);
		double[] numbers = new double[strings.length];
		for (int i = 0; i < strings.length; i++)
		{
			numbers[i] = Double.parseDouble(strings[i]);
		}
		return numbers;
	}

	private static double sigmoid(double value)
	{
		return (1.0 / (1.0 + Math.exp(-value)));
	}
}